Numerically generated quasi-equilibrium orbits of black holes: Circular or eccentric? 
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We make a comparison between results from numerically generated, quasi-equilibrium configura- 
tions of compact binary systems of black holes in close orbits, and results from the post-Newtonian 
approximation. The post-Newtonian results are accurate through third PN order (0(v/c) 6 beyond 
Newtonian gravity), and include rotational and spin-orbit effects, but are generalized to permit 
orbits of non-zero eccentricity. Both treatments ignore gravitational radiation reaction. The energy 
E and angular momentum J of a given configuration are compared between the two methods as 
a function of the orbital angular frequency fl. For small Q, corresponding to orbital separations a 
factor of two larger than that of the innermost stable orbit, we find that, if the orbit is permitted 
to be slightly eccentric, with e ranging from ~ 0.03 to » 0.05, and with the two objects initially 
located at the orbital apocenter (maximum separation), our PN formulae give much better fits to 
the numerically generated data than do any circular-orbit PN methods, including various "effective 
one-body" resummation techniques. We speculate that the approximations made in solving the 
initial value equations of general relativity numerically may introduce a spurious eccentricity into 
the orbits. 



I. INTRODUCTION AND SUMMARY 

The late stage of inspiral of binary systems of neu- 
tron stars or black holes is of great current interest, both 
as a challenge for numerical relativity, and as a possible 
source of gravitational waves detectable by laser inter- 
feromctric antennas. Because this stage, corresponding 
to the final few orbits and ultimate merger of the two 
objects into one, is highly dynamical and involves strong 
gravitational fields, it must be handled by numerical rela- 
tivity, which attempts to solve the full Einstein equations 
on computers (see Refs. [fil || for reviews). 

The early stage of inspiral can be handled accurately 
using post-Newtonian techniques, which involve an ex- 
pansion of solutions of Einstein's equations in powers of 
e ~ (v/c) 2 ~ Gm/rc 2 , where v, m, and r are typical ve- 
locity, mass and separation in the system, respectively. 
By expanding to very high powers of e, increasingly ac- 
curate formulae can be derived to describe both the or- 
bital motion and the gravitational waveform. Currently, 
results accurate through 3.5PN order (CHe 7 ^ 2 ) beyond 
Newtonian gravity) are known || f§ §, g, @, |, g. 

An important issue in understanding the full inspiral of 
compact binaries is how to connect the PN regime to the 
numerical regime. This is a non-trivial issue because the 
PN approximation gets worse the smaller the separation 
between the bodies. On the other hand, because of lim- 



tmora@cliDDcr.cns.fr 


:mw@wuphys. wustl.edu 


»ugr av . wustl . edu/People/CLIFF 



ited computational resources, numerical simulations can- 
not always be started with separations sufficiently large 
to overlap the PN regime where it is believed to be reli- 
able. 

Numerical simulations of compact binary inspiral start 
with a solution of the initial value (I-V) equations of Ein- 
stein's theory; these provide the initial data for the evo- 
lution equations (some simulations ]To| solve in addition 
one of the six dynamical field equations) . The initial state 
is assumed to consist of two compact objects (neutron 
stars or black holes) in an initially circular orbit. The 
circular-orbit condition is imposed by demanding that 
dr/dt = initially, where r is a measure of the orbital 
separation. More precisely, the system is presumed to 
have an initial "helical Killing vector" (HKV), which cor- 
responds to a kind of rigid rotation of the binary system. 
This amounts to ignoring initially the effects of gravita- 
tional radiation reaction. It also implies that the black 
holes are co-rotating, a condition which is astrophysically 
unlikely, albeit computationally advantageous. To make 
the problem tractable numerically, it is also generally as- 
sumed that the spatial metric is conformally flat. This 
approximation is usually justified by the neglect of radia- 
tion reaction in the initial state. For black hole binaries, 
suitable horizon boundary conditions must be imposed, 
while for neutron star binaries, hydrodynamical equa- 
tions and an equation of state must be provided. 

One important product of these initial value solutions 
is a relationship between the energy E and angular mo- 
mentum J of the system as measured at infinity, and the 
orbital frequency f2. As all quantities are well-defined 
and gauge invariant, they are useful variables for making 
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comparisons with PN methods. 

We have developed a formula for E(d) and J(0) using 
PN methods. Our analytic formula includes point-mass 
terms through 3PN order, but ignores radiation reaction. 
It also includes rotational energy and spin-orbit terms for 
the case in which the bodies are rotating. For black holes, 
tidal effects can be ignored. In contrast to previous work 
p[ |l2[ , our formulae apply to general eccentric orbits, 
not just to circular orbits. 

We then compare this formula with HKV numerical 
solutions for corotating binary black holes obtained by 
Grandclement et al. [Q, for the regime where the black 
holes are separated from the location of the innermost 
circular orbit by a factor of around two, where PN results 
might be expected to work well (Gm/rc 2 ~ 0.1). We find 
the following two results: 

1. When we assume circular PN orbits, our PN formulae 
for E(£l) and J(r2) agree well with other PN methods, 
including those using resummation or Pade techniques. 
However all PN methods consistently and systematically 
underestimate the binding energy and overestimate the 
angular momentum, compared to the values derived from 
the numerical HKV simulations, by amounts that are 
larger than the spread among the PN methods. 

2. When we relax the assumption of a circular orbit and 
demand only that dr/dt = 0, our PN formula agrees ex- 
tremely well with the numerical data. In this case the sys- 
tem is found to be initially at the apocenter of a slightly 
eccentric orbit. For values of Gmfl/c 3 ranging from 0.03 
to 0.06, corresponding to orbital v/c between 0.3 and 0.4, 
or orbital separation between 10 and 6 Gm/c 2 , perfect 
agreement with the binding energy can be obtained with 
eccentricities that range from 0.03 to 0.05. 

If this apparent eccentricity is a real effect, there are 
two possible explanations. One is that the imposition of 
a helical Killing vector, which is equivalent to assuming a 
stationary binary star configuration as seen in a rotating 
frame, is not sufficient to guarantee a circular orbit in the 
relativistic regime. In Newtonian gravity, this assump- 
tion would be equivalent to imposing dr/dt = 0, which 
only puts the orbit at a turning point. An additional con- 
dition d 2 r/dt 2 = is needed to enforce a circular orbit. 
But this is a dynamical statement, which is outside the 
realm of the I-V equations of Einstein's theory. It may 
be that, in solving one of the dynamical field equations 
(the trace of the equation for dK 1 ^ /dt) in addition to the 
four I-V equations, Grandclement et al. take care of this 
automatically. 

The more likely possibility is that the approximations 
made in most numerical solutions, the main one being 
that the spatial metric is conformally flat, somehow in- 
troduce a systematic eccentricity into the orbit. 

At present, these results can only be considered a hint 
of a possible eccentricity, however. It is entirely possible 
that circular PN orbits are consistent with the HKV re- 
sults within the errors of the numerical simulations. This 
can only be decided if and when the numerical groups 
that carry out these simulations publish quantitative er- 



ror bars determined by studying the convergence prop- 
erties of the solutions as a function of grid size, domain 
size and other computational assumptions. 

The remainder of this paper sketches the arguments 
that lead to these conclusions. Detailed derivations and 
formulae, and applications to neutron-star simulations, 
will be the subject of a future publication. 



II. ORBITS AT THE TURNING POINT IN 
POST-NEWTONIAN GRAVITY 

In Newtonian gravity, the orbit of a pair of point 
masses may be described by the set of equations 
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1 + ecos(0 — uo) , 
r 2 d(f>/dt = (mp) 1/2 , 
fi(r 2 +r 2 tt 2 )/2- fim/r, 
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where p — a(l — e 2 ) is the semi-latus rectum (a is the 
semi- major axis), uo is the angle of pericenter, m = 
ui\ + rri2 is the total mass, fi = m\mijm is the re- 
duced mass, and E is the total orbital energy (hence- 
forth we use units in which G = c = 1). A circular 
orbit corresponds to e = 0, with r = a = constant, 
n 2 = m/a 3 , E/n = a 2 n 2 /2 - m/a = -(mQ) 2 / 3 /2, and 
J / \x = y/ma — (m/fi) 1 / 3 . However, if we demand only 
that the orbit be at apocenter, so that r = only, we 
have <p = uo + tt, r = p/(l — e), ft 2 = (m/p 3 )(l — e) 4 , and 
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where Q a is the angular velocity at apocenter. The en- 
ergy and angular momentum expressions when the orbit 
is at pericenter may be obtained by letting f2 a 



rip and 



At 1PN order, the equations of motion can be solved 
using either the direct approach of Wagoner and Will 
]l3[ , or the "perturbed osculating orbit elements ap- 
proach" of Lincoln and Will |l4| to yield the orbit equa- 
tions 
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1 + e cos(</> — uj) 

+C[-(3-77) + (l + 9?7/4)g 2 

+ (7/2 — ri)ecos(4> — uj) + 3e0sin(0 — uj) 

-(?7/4)e 2 cos(20 - 2uj)] + O(C) 2 , 

(mp) 1/2 [l - C(4 - 2r 1 )icos((t) - uj) + 0(C) 2 ] 

fi(r 2 +r 2 D, 2 )/2 - fim/r 

3 11 

+ -fj,v 4 (l - 3r?) + -[i(m/r) 2 + -^rjr 2 m/r 
o 2 2 

3 

+ -[i(v 2 m/r)(l + 77/3) + 0(m/r) 4 , 
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where v 2 = r 2 + r 2 fi 2 , r\ = /i/m, and £ = m/p. The limit 
e — > corresponds to a circular PN orbit. However, at 
higher PN orders, neither the orbital eccentricity e nor 
the semilatus rectum p is uniquely or invariantly defined. 
One definition of eccentricity used by Lincoln and Will 
fhij was that of a Newtonian orbit momentarily tangent 
to the true orbit (the "osculating" eccentricity); other au- 
thors [[l5| define multiple "eccentricities" to encapsulate 
different aspects of non-circular orbits at PN order. 

We define an alternative eccentricity and semilatus rec- 
tum according to: 
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where f2 p is the value of f2 where it passes through a local 
maximum (pericenter), and £l a is the value of f2 where 
it passes through the next local minimum (apocenter). 
These definitions have the following virtues: (1) they re- 
duce precisely to the normal eccentricity e and semilatus 
rectum p in the Newtonian limit, as can be seen from Eqs. 
(|l|); (2) they are constant in the absence of radiation re- 
action; (3) they are somewhat more directly connected to 
measurable quantities, since is the angular velocity as 
seen from infinity (eg. as measured in the gravitational- 
wave signal) and one calculates only maximum and mini- 
mum values, without concern for the coordinate location 
in the orbit; and (4) they are straightforward to calcu- 
late in a numerical simulation of orbits without resorting 
to complicated definitions of "distance" between bodies. 
They have the defect that, when radiation reaction is 
included, they are not local, continuously evolving vari- 
ables, but rather are some kind of orbit-averaged quan- 
tities (for this reason, they may not be as "covariant" 
as they seem - this issue is under investigation). Never- 
theless, when an eccentric orbit decays and circularizes 
under radiation reaction the definition of e has the virtue 
that it tends naturally to zero when the orbital frequency 
turns from ocillatory behavior to monotonically increas- 
ing behavior (i.e. the maxima and minima merge). 

By virtue of these definitions, £ has the further prop- 
erty that 
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It is then simple to show that the relation between e and 
£ and the corresponding quantities used in the Wagoner- 
Will solution (§ is e = e{\ + C[9/2 - rj + (1 - 3r/)e 2 ]}, 
and C = C{1 - 4C[1 - 7/73 + (1/3 - v)i 2 }}. Applying 
these definitions to Eqs. (^), we find the 1PN expression 



for E(fl) and J(fi) for a non-circular orbit, expressed in 
terms of f2 at the apocenter or pericenter: 
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Note that in the circular orbit limit (e — > 0), E and 
J satisfy the expected relation dE/dfl = fldJ/dil. We 
have extended this result to 3PN order using the 3PN 
orbit equations of Blanchet and Faye @, using both 
harmonic-gauge and ADM gauge expressions, but ne- 
glecting radiation-reaction terms at 2.5PN order (details 
will be given in a future publication). 

For those numerical simulations in which the bodies 
are assumed to be corotating, we must include a num- 
ber of rotational effects. First is the energy of rotation 
of the individual corotating bodies, of order SE TOt /p, ~ 
Ift 2 /fi ~ (m/r){R/r) 2 , where R is the characteristic size 
of the body. For compact bodies, R ~ M, so these ef- 
fects are equivalent to 2PN order and higher. For black 
holes, we will use the standard procedure of splitting the 
mass of each body into its irreducible mass and its ro- 
tational energy using the Kerr formulae M = Mj rr /[1 — 
4(M irr fi) 2 ] 1 / 2 and S = 4M3. r ft/[l-4(M rr ft) 2 ] 1 / 2 00- 
Also, since the sequences of numerical simulations hold 
m m = (Mj rr )i + (Mi rr )2 fixed as they vary f2, we will 
expand the masses that appear in the Newtonian orbital 
contribution to the energy and angular momentum about 
rrii rr . These contributions together yield 

SErat/lt = (2/r/)(l-377)(m irr f>) 2 

+ (6/77)(l- 577 + 577 2 )(m irr ft) 4 

-(1 - e 2 )(l - 77)C(m irr f!) 2 , 



SJrot/lMTli, 



(4/77)(l-3r/)(m lrr O) 

+ (8/7 ? )(l-5r ? + 5r? 2 )(m irr O) 3 



1 ~N2,4 



+ -j={m iTT ny(~-2r)).. 



(7) 



where henceforth, £ = [mi rr f2 a /(l — e) 2 ] 2 / 3 , and fi and 
77 are expressed in terms of irreducible masses. The first 
two terms in each of Eqs. (^) are the rotational terms, 
expanded to second order, while the third comes from 
expanding the Newtonian orbital term. 

Similarly spin-orbit and spin-spin effects are of order 
SEs.o./fJ- ~ LSj [ir z and SEs.s./fJ- ~ SiS 2 / p.r 3 where L 
denotes the orbital angular momentum and S denotes 
the bodies' spin. For co-rotating bodies, these effects 
can be shown to be of order SEg.o./l- 1 ~ (m/r)(mR 2 /r 3 ) 
and SEs.s./n ~ (m/r)(mR 4 /r 5 ), which for compact bod- 
ies are equivalent to 3PN and 5PN order, respectively. 
Henceforth, we will ignore the spin-spin terms. Assum- 
ing co-rotating bodies with spins aligned with the orbital 
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FIG. 1: Binding energy of equal-mass, co-rotating black holes 
vs. angular frequency. 



FIG. 2: Angular momentum of equal-mass, co-rotating black 
holes vs. angular frequency. 



angular momentum, including both the direct and or- 
bital effects of the spin-orbit terms jl6|, [lTj as well as 
their effect on our definitions of e and £, we find 
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Apart from the generalization to eccentric orbits, these 
methods parallel exactly those of Blanchet jll] . 

For tidal interactions, the contribution to the orbital 
energy is of order #-E?tidai//-t ~ (m/r)(R/r) 5 . For black 
holes, tidal effects are thus equivalent to 5PN order, and 
will be neglected. For neutron stars, R ~ 5M, and tidal 
effects must be included, however it is sufficient to use 
Newtonian gravity to calculate them. These will be dis- 
cussed elsewhere. 
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FIG. 3: Values of eccentricity giving a match to numerical 
HKV simulations. 



III. COMPARISON WITH NUMERICAL 
SIMULATIONS OF COROTATING BINARY 
BLACK HOLES 

We combine the PN, rotational and spin-orbit contri- 
butions to E and J, set fi a = so that the stars corotate 
at the orbital angular frequency at apocenter, set 77 = 1/4 
for equal-mass black holes, and plot (E — mi 1T )/mi rr and 
J/mj rr versus mi rr f2, where mi rr denotes the total irre- 
ducible mass. The results are shown in Figures 1 and 
2. For e = 0, we show results both from the 3PN or- 
bital expressions, and from truncated 2PN orbital ex- 
pressions, while for e = 0.03 we plot the full 3PN results. 
The crosses denote the corresponding numerical data of 
Grandclement et al.p0[; the other curve is the 3PN re- 
summation result of Damour et al.[fl2"| for circular orbits. 
We note that all circular-orbit PN estimates, including 



conventional 2PN and 3PN approaches JO], Pade and 
resummation approaches |p^| , and our approach, agree 
well among themselves, but are systematically displaced 
from the numerical points, both for E and for J. This 
can also be seen in Figures 5 and 6 of |12] |, By contrast, 
with e ~ 0.03, our PN estimate fits the curves perfectly 
for small values of m m VL. This fit is merely illustrative: 
if the numerical simulations do have some eccentricity, 
there is no reason why the eccentricity should be the 
same for each numerical point. Figure 3 shows the val- 
ues of e that best fit the numerical data over the range 
of Q, considered. 

The conclusions are not substantially different if we 
use the Wagoner- Will type eccentricity e. Whatever def- 
inition one uses, e = gives a worse fit than a finite 
eccentricity. 
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FIG. 4: Comparison of PN calculations of energy for widely 
separated black holes. 



IV. DISCUSSION 

High-order post-Newtonian approximations for circu- 
lar orbits appear to give results for the energy and angu- 
lar momentum of corotating binary black holes well away 
from the innermost orbit that are in excellent agreement 
- with each other. This is illustrated in Figure 4, which 
plots the energy from four PN results quoted by Damour 
et al JT^ ], and from our 2PN and 3PN results, in the 
range around Q,m- m ~ 0.32. Apart from the one 2PN re- 
sult quoted by Damour et al., all are within 0.5 percent 
of each other. This could be viewed as an estimate of the 
accuracy of the PN approximation in this regime. But 
all results are four percent displaced from the numerical 



HKV data. In the absence of a quantitative estimate of 
the accuracy of the numerical simulations in this regime, 
it is difficult to decide if this difference is a signal of a 
physical effect, such as the small eccentricity suggested 
by our work. 

One way to test whether the orbits represented by the 
HKV numerical simulations are really eccentric would 
be to evolve the orbits numerically for a short period of 
time, say 1/4 of an orbit (5(f) = it/2). For a nearly cir- 
cular orbit of equal mass bodies, gravitational radiation 
reaction should cause the parameter £ — m/p = (mfl) 2 / 3 
to evolve according to dC,/d(j) = (16/5)C 7/2 (see, for exam- 
ple Eq. (3.11a) of 0). Thus, over a quarter of an orbit, 
C should change by approximately <5C/C ~ (87r/5)£ 5 ' 2 « 
(87r/5)(m^) 5 / 3 w 0.016, for mil ~ 0.032. Hence the or- 
bital separation should decrease by about 1.6 percent or 
the angular velocity should increase by 2.4 percent. By 
contrast, if the orbit has an eccentricity of 0.03 and is at 
apocenter, then in a quarter of an orbit, the separation 
should decrease by 3 percent, while the angular velocity 
should increase by 6 percent. Even with radiation damp- 
ing, the orbit should pass through a distinct pericenter 
when 84> ~ 7r, and the orbital separation should increase, 
while the angular velocity decreases. 
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